library(dplyr)
library(tidyr)
library(ggplot2)
library(lubridate)
library(data.table)
library(purrr)
library(stringr)
library(plm)
library(lmtest)
library(sandwich)

##############
## Analysis ##
##############
pdata <- read.csv("acled_local_final.csv")

library(slider)
library(lubridate)

pdata <- pdata %>%
  mutate(event_date = as.Date(event_date)) %>%
  arrange(admin2, event_date) %>%
  group_by(admin2) %>%
  mutate(
    n_events_recent365 = slide_index_int(
      .x = event_date,    
      .i = event_date,    
      .f = ~length(.x),   
      .before = days(365),  
      .complete = FALSE
    )
  ) %>%
  ungroup()


pdata <- pdata.frame(pdata, index = c("admin2", "event_date"))

pdata_quad <- pdata %>%
  filter(
    !is.na(fatality_diff_signed_weighted),
    !is.na(participants_lead),
    !is.na(fatalities_log),
    !is.na(protest),
    !is.na(civilian)
  )

poly_terms <- poly(pdata_quad$fatality_diff_signed_weighted, 2)

m_quad <- plm(
  participants_lead ~ poly(fatality_diff_signed_weighted, 2) + n_events_recent365 +
      protest + civilian + v2x_polyarchy + v2xel_locelec,
  data = pdata_quad,
  index = c("admin2", "event_date"),
  model = "within",
  effect = "twoways"
)

acled_result1_admin <- coeftest(m_quad, vcov = vcovHC(m_quad, type = "HC1", cluster = "group"))

fatality_range <- seq(
  min(pdata_quad$fatality_diff_signed_weighted, na.rm = TRUE),
  max(pdata_quad$fatality_diff_signed_weighted, na.rm = TRUE),
  length.out = 100
)

poly_deriv <- predict(poly_terms, newdata = fatality_range, deriv = 1)

coefs <- coef(m_quad)
vcov_mat <- vcovHC(m_quad, type = "HC1", cluster = "group")
results3 <- coeftest(m_quad, vcov. = vcov_mat)

marginal_effect <- poly_deriv[, 1] * coefs["poly(fatality_diff_signed_weighted, 2)1"] +
  poly_deriv[, 2] * coefs["poly(fatality_diff_signed_weighted, 2)2"]

grad_matrix <- cbind(poly_deriv[, 1], poly_deriv[, 2])
colnames(grad_matrix) <- c("poly(fatality_diff_signed_weighted, 2)1",
                           "poly(fatality_diff_signed_weighted, 2)2")

marginal_se <- apply(grad_matrix, 1, function(g) {
  sqrt(t(g) %*% vcov_mat[names(g), names(g)] %*% g)
})

prediction_data <- data.frame(
  fatality_diff_signed_weighted = fatality_range,
  marginal_effect = marginal_effect,
  se = marginal_se
) %>%
  mutate(
    ci_low = marginal_effect - 1.96 * se,
    ci_high = marginal_effect + 1.96 * se
  )

acled_marginal1_admin <- ggplot(prediction_data, aes(
  x = fatality_diff_signed_weighted, y = marginal_effect
)) +
  geom_line(color = "black", size = 0.5) +
  geom_ribbon(aes(ymin = ci_low, ymax = ci_high), fill = "black", alpha = 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Marginal Effect of Repression Deviations\n(ACLED: Local)",
    x = "Repression Deviations",
    y = "Marginal Effect on Protest Participation (log)"
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 14),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    axis.title.x = element_text(size = 15, margin = margin(t = 15)),
    axis.title.y = element_text(size = 15, margin = margin(r = 15)),
    axis.text = element_text(size = 15)
  )

print(acled_marginal1_admin)

pdata_quad <- pdata %>%
  filter(
    !is.na(fatalities_log),
    !is.na(participants_lead),
    !is.na(protest),
    !is.na(v2xel_locelec),
    !is.na(v2x_polyarchy)
  )

poly_terms <- poly(pdata_quad$fatalities_log, 2)

m_quad_alt <- plm(
  participants_lead ~ poly(fatalities_log, 2) + protest + civilian + n_events_recent365 +
  + v2x_polyarchy + v2xel_locelec,
  data = pdata_quad,
  index = c("admin2", "event_date"),
  model = "within",
  effect = "twoways"
)

acled_result2_local <- coeftest(m_quad_alt, vcov = vcovHC(m_quad_alt, type = "HC1", cluster = "group"))

fat_range <- seq(
  min(pdata_quad$fatalities_log, na.rm = TRUE),
  max(pdata_quad$fatalities_log, na.rm = TRUE),
  length.out = 100
)

poly_deriv <- predict(poly_terms, newdata = fat_range, deriv = 1)

coefs <- coef(m_quad_alt)
vcov_mat <- vcovHC(m_quad_alt, type = "HC1", cluster = "group")

marginal_effect <- poly_deriv[, 1] * coefs["poly(fatalities_log, 2)1"] +
  poly_deriv[, 2] * coefs["poly(fatalities_log, 2)2"]

grad_matrix <- cbind(poly_deriv[, 1], poly_deriv[, 2])
colnames(grad_matrix) <- c("poly(fatalities_log, 2)1", "poly(fatalities_log, 2)2")

marginal_se <- apply(grad_matrix, 1, function(g) {
  sqrt(t(g) %*% vcov_mat[names(g), names(g)] %*% g)
})

prediction_data <- data.frame(
  fatalities_log = fat_range,
  marginal_effect = marginal_effect,
  se = marginal_se
) %>%
  mutate(
    ci_low = marginal_effect - 1.96 * se,
    ci_high = marginal_effect + 1.96 * se
  )

acled_marginal2_admin <- ggplot(prediction_data, aes(x = fatalities_log, y = marginal_effect)) +
  geom_line(color = "black", size = 0.5) +
  geom_ribbon(aes(ymin = ci_low, ymax = ci_high), fill = "black", alpha = 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Marginal Effect of Simple Repression\n(ACLED: Local)",
    x = "Fatalities",
    y = "Marginal Effect on Protest Participation (log)"
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 14),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    axis.title.x = element_text(size = 15, margin = margin(t = 15)),
    axis.title.y = element_text(size = 15, margin = margin(r = 15)),
    axis.text = element_text(size = 15)
  )

print(acled_marginal2_admin)

# Step 1: Filter data
pdata_quad_exp <- pdata %>%
  filter(
    !is.na(fatality_diff_signed_weighted_exp),
    !is.na(participants_lead),
    !is.na(fatalities_log),
    !is.na(protest),
    !is.na(civilian)
  )

# Step 2: Fit quadratic plm model with exponential variable
poly_terms_exp <- poly(pdata_quad_exp$fatality_diff_signed_weighted_exp, 2)

m_quad_exp <- plm(
  participants_lead ~ poly(fatality_diff_signed_weighted_exp, 2) + n_events_recent365 +
    protest + civilian + v2x_polyarchy + v2xel_locelec,
  data = pdata_quad_exp,
  index = c("admin2", "event_date"),
  model = "within",
  effect = "twoways"
)

acled_result1_admin_exp <- coeftest(m_quad_exp, vcov = vcovHC(m_quad_exp, type = "HC1", cluster = "group"))

# Step 3: Create range of repression deviation values
fatality_range_exp <- seq(
  min(pdata_quad_exp$fatality_diff_signed_weighted_exp, na.rm = TRUE),
  max(pdata_quad_exp$fatality_diff_signed_weighted_exp, na.rm = TRUE),
  length.out = 100
)

poly_deriv_exp <- predict(poly_terms_exp, newdata = fatality_range_exp, deriv = 1)

coefs_exp <- coef(m_quad_exp)
vcov_mat_exp <- vcovHC(m_quad_exp, type = "HC1", cluster = "group")

marginal_effect_exp <- poly_deriv_exp[, 1] * coefs_exp["poly(fatality_diff_signed_weighted_exp, 2)1"] +
  poly_deriv_exp[, 2] * coefs_exp["poly(fatality_diff_signed_weighted_exp, 2)2"]

grad_matrix_exp <- cbind(poly_deriv_exp[, 1], poly_deriv_exp[, 2])
colnames(grad_matrix_exp) <- c("poly(fatality_diff_signed_weighted_exp, 2)1",
                               "poly(fatality_diff_signed_weighted_exp, 2)2")

marginal_se_exp <- apply(grad_matrix_exp, 1, function(g) {
  sqrt(t(g) %*% vcov_mat_exp[names(g), names(g)] %*% g)
})

# Step 4: Prepare data for plot
prediction_data_exp <- data.frame(
  fatality_diff_signed_weighted_exp = fatality_range_exp,
  marginal_effect = marginal_effect_exp,
  se = marginal_se_exp
) %>%
  mutate(
    ci_low = marginal_effect - 1.96 * se,
    ci_high = marginal_effect + 1.96 * se
  )

# Step 5: Plot
acled_marginal1_admin_exp <- ggplot(prediction_data_exp, aes(
  x = fatality_diff_signed_weighted_exp, y = marginal_effect
)) +
  geom_line(color = "black", size = 0.5) +
  geom_ribbon(aes(ymin = ci_low, ymax = ci_high), fill = "black", alpha = 0.3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  labs(
    title = "Marginal Effect of Repression Deviations\n(ACLED: Country, Exponential Decay)",
    x = "Repression Deviations (Exponential Decay)",
    y = "Marginal Effect on Protest Participation (log)"
  ) +
  theme_minimal() +
  theme(
    text = element_text(size = 14),
    plot.title = element_text(face = "bold", size = 16, hjust = 0.5),
    axis.title.x = element_text(size = 15, margin = margin(t = 15)),
    axis.title.y = element_text(size = 15, margin = margin(r = 15)),
    axis.text = element_text(size = 15)
  )

# Step 6: Print plot
print(acled_marginal1_admin_exp)

